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ABSTRACT 

Accretion onto black holes and compact stars brings material in a zone of strong gravitational 
and electromagnetic fields. We study dynamical properties of motion of electrically charged particles 
forming a highly diluted medium (a corona) in the regime of strong gravity and large-scale (ordered) 
magnetic field. 

We start our work from a system that allows regular motion, then we focus on the onset of chaos. To 
this end, we investigate the case of a rotating black hole immersed in a weak, asymptotically uniform 
magnetic field. We also consider a magnetic star, approximated by the Schwarzschild metric and a test 
magnetic field of a rotating dipole. These are two model examples of systems permitting energetically 
bound, off-equatorial motion of matter confined to the halo lobes that encircle the central body. Our 
approach allows us to address the question of whether the spin parameter of the black hole plays any 
major role in determining the degree of the chaoticness. 

To characterize the motion, we construct the Recurrence Plots (RP) and we compare them with 
Poincare surfaces of section. We describe the Recurrence Plots in terms of the Recurrence Quantifica- 
tion Analysis (RQA), which allows us to identify the transition between different dynamical regimes. 
We demonstrate that this new technique is able to detect the chaos onset very efficiently, and to 
provide its quantitative measure. The chaos typically occurs when the conserved energy is raised to a 
sufficiently high level that allows the particles to traverse the equatorial plane. We find that the role 
of the black- hole spin in setting the chaos is more complicated than initially thought. 
Subject headings: acceleration of particles, black hole physics, magnetic fields, methods: numerical 



1. INTRODUCTION 

The role of magnetic fields near strongly gravitating 
objects has been subject of many investigations (e.g. 
iPunslvl [2008). They are relevant for accretion disks 
that may be embedded in large-scale magnetic fields, 
for example whe n the accretion flow penetrat es close to 
a neutron star ([Lipunovl 119921 : iMestell fl999f ). Outside 
the main body of the accretion disk, i.e. above and be- 
low the equatorial plane, the accreted material forms a 
highly diluted environment, a 'corona', where the den- 
sity of matter is low and the mean free path of parti- 
cles is large in comparison with the characteristic length- 
scale, i.e. the gravitational radius of the central body, 
R g GM/c 2 « 1.5(M/M ) km, where M is the cen- 
tral mass. The origin of the coronal flows and the rele- 
vant processes governing their structure are still unclear. 
In this context we discuss motion of electrically charged 
particles outside the equatorial plane. 

Off-equatorial, energetically bound motion of charged 
particles in strong gravitational and electromagnetic 
fields is pertinent to the description of accretion disk 
coronae around black h oles and compact (neutro n) stars. 
In our previous papers (jKovaf et al.l 12008. 2010) we dis- 
cussed the existence of energetically-bound stable or- 
bits of charged particles occurring outside the equato- 
rial plane, extending thus a large variety of comple- 
mentary studies (e.g.. [Bleak et al][l98 9: Stuc hlik et al.l 



119991: [Prasannal 119801 : iPrasanna fe Senguptal 119941 : 
lAliev fc Ozdemirl I2001L and further references cited 
therein). Particles on off-equatorial stable trajectories 
form a coronal flow that is possible at certain radii 
and for certain combinations of the model parameters, 
namely, the specific charge of the particles, the conserved 
energy and the angular momentum of the particle mo- 
tion, the strength and orientation of the magnetic field, 
and the spin of the central body. 

We assume that the magnetic field permeating 
the corona has a large-sca l e (or dered) component 
(Bis novatvi-Kogan fc Lovelacd 120071 ). In this case, 
charged particles can be trapped in toroidal regions, ex- 
tending symmetrically above and below the equatorial 
plane and forming two halo lobes. However, this trap- 
ping happ ens only for certai n combinations of model pa- 
rameters (|Kovaf et al.ll2010T ). 

We consider two types of the model setup: a ro- 
tating (Kerr) black hole in an asymptotically uni- 
form magnetic field parallel to the symmetry axis 
(lWaldlll974t iTomimatsu fc Takahashil 120011 : [Koidd 120041 : 
iKoide et al. 2006), and a non-rotating star (described by 
the Schwarzschild metric) endowed with a r otating mag- 
netic dipole field (iPettersonll 19751 : lSenguptalll995l ). Both 
cases can be regarded as integrable systems with the elec- 
tromagnetic field acting as a perturbation. 

The above-mentioned lobes are defined by the figures 
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Figure 1. In the left panel we present a poloidal section of the selected isocontours of the effective potential V e ^(r, 9), eq. (0), for a charged 
particle (qQ = 2, L = 5 M) on the Kerr background (a = 0.5 M). We assume the presence of Wald uniform magnetic field (qBo = 2M _1 ). 
The off-equatorial potential lobes are present, allowing stable motion. Two exemplary trajectories of test particles are shown - in the left 
lobe a chaotic orbit of energy E = 1.796, while in the right lobe the regular, purely off-equatorial trajectory of E = 1.78. Both particles 
were launched at r(0) = 3.11, 0(0) = 7r/4 with u r (0) = and their trajectories interweave with each other. We plot the poloidal (r, 9) 
projection of the trajectory; what appears as a lobe in the poloidal plane is an axially symmetric 3-dimensional rotational structure. The 
latter is illustrated in the right panel where the case of the off-equatorial regular trajectory is shown. 



of the effective potential in the poloidal plane. These 
were previously studied in the context of charge separa- 
tion that is expecte d to occur in pulsar magnetospheres 
(e.g. lNeukirch|[l993l ). Here, we address whether the tra- 
jectories within these lobes are regular (i.e., whether the 
system is integrable), or if they instead exhibit a chaotic 
behavior. A related proble m was studied recently by 
Takahash ifc Kovamal (2009) in an attempt to find a con- 
nection between chaoticness of the motion and the spin 
of a rotating black hole residing in the center. These 
authors suggest that chaotic behavior occurs for certain 
values of the black hole spin, while for others the system 
is indeed regular. 

The idea of investigating the connection between the 
spin of a black hole and chaoticness of motion of mat- 
ter near its horizon is very interesting for the follow- 
ing reason. Because of high degree of symmetry of the 
back ground spac etime, the unperturbed motion is regu- 
lar ([Carter! fl968): no chaos is present. The electromag- 
netic perturbation may trigger the chaos, however, its 
effect can be expected to diminish very near the horizon, 
where strong gravity of the black hole should prevail. 
This is also the region where the spin effects are most 
prominent. Further out various other influences become 
important due to distant matter and the turbulence in 
accreted material. Therefore the connection between the 
spin and the motion chaoticness is best applicable in the 
immediate vicinity of the black hole, i.e. within the inner 
parts of corona. 

The recurrence analysis (M arwan et al.l l20Q7) provides 
us with a powerful tool for the investigation of complex 
dynamical systems. The method examines the recur- 
rences of the system to the vicinity of previously reached 
phase space points. It has been typically adopted to 
study the experimental data, where often only some (if 
not just one) of the phase space variables are kno wn from 
the m easurements. Takens' embedding theorems (jTakens! 
Il981[ ) are then used to reconstruct the phase space por- 



trait of such a system. In our study we are equipped 
with the full phase space trajectory from the numerical 
integration of the equations of motion, so that we can 
use the recurrence analysis directly. 

It appears that the method of Recurrence Plots has not 
been employed in the context of relativistic astrophysical 
systems yet. To this end, one needs a consistent defini- 
tion of the neighborhood of a point in the phase space in 
a curved spacetime. Below, we discuss the phase space 
distance and suggest a form of the distance norm suitable 
in such circumstances. 

The paper is organized as follows. In sec. [2] we review 
the equations of particle motion, which we then integrate 
to obtain trajectories. In sec. [3] we introduce the basic 
properties of Recurrence Plots. Sec. [4] analyses the mo- 
tion around a Kerr black hole endowed with a uniform 
magnetic test field. We employ Poincare surfaces of sec- 
tion and Recurrence Plots. The two approaches allow us 
to show the onset of chaos in different, complementary 
ways. We examine the motion in off-equatorial lobes, pay 
special attention to the spin dependence of the stability 
of motion, and we notice the emergence of 'potential val- 
leys' that allow the particles to escape from the equato- 
rial plane along a narrow collimated corridor. Analysis of 
the off-equatorial motion around a magnetic star is pre- 
sented in sec. [51 We consider a dipole-type magnetic field, 
which sets different limits on the off-equatorial range of 
allowed motion of charged particles. It also defines dif- 
ferent regimes of chaoticness in comparison with the uni- 
form magnetic field. Finally, results of the analysis are 
summarized in sec. [U 

2. EQUATIONS OF MOTION AND THE EFFECTIVE 
POTENTIAL 

The phase space trajectories of integrable systems are 
regular, meaning that they are bound to the surface of 
an n-dimensional torus, where n is the number of de- 
grees of freedom. The torus is determined uniquely by 
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Figure 2. The overview of possible topologies of the off-equatorial potential structure above the event horizon (thick line in plots) of Kerr 
black hole endowed with the Wald test field. 



n constants of motion that are present in such a sys- 
tem. Its behavior can be explored by Poincare sur- 
faces of section, which are defined by intersections of 
the phase space trajectory with a 2-dimensional plane 
(jLichtenberg fc Liebermanl Il992l ). On the other hand, 
non-integrable chaotic systems generally have fewer inte- 
grals than the number of degrees of freedom. In general, 
both the regular and the chaotic orbits may coexist in 
the phase space of a single system. 

Chaotic orbits are ergodic on the given hypersurface. 
Its dimension is now larger than n, and the section points 
thus fill areas in the plot of the Poincare surface. How- 
ever, depending on the initial conditions, regular orbits 
can also appear in non-integrable systems. Such orbits 
maintain the value of some additional constant of motion, 
although it is not generally possible to write this constant 
in an explicit form. In the context of motion around black 
holes perturbed by (weak) exter nal sources, various as- 
pects o f chaos were studied e.g . bvlkaras fc Vok rouhlickv 
(1992); Naka mura fc Ishizu ka (1993); Podolskv fc Veselvl 
(1998), and very recently bv lSemerak fc Sukoval (|2010l ). 

A standard approach to an integrable system with 
a non-integrable perturbation assumes complete control 
over the strength of the perturbation (i.e., the pertur- 
bation can be set to be arbitrarily weak). If this were 
the case, we could first switch the perturbation com- 
pletely off, analyze the orbits, and then observe the im- 
pact of gradually increasing the perturbation strength 
upon these orbits. However, the class of off-equatorial 
bound orbits only exists when the electromagnetic term 
is strong enough to balance the gravitational attraction 
of the central body. Then the (sufficiently strong) pertur- 
bation is by itself the cause of the new kind of the regular 
motion that happens outside the equatorial plane. 

Having this delicacy on mind, we shall use the usual 
Hamiltonian formalism to express equations of motion 
governing the tr ajectories. We first construct the super- 
Hamiltonian U (|Misner et al.lll973D . 

(1) 



g^ v is the metric tensor, and denotes the vector po- 
tential of the electromagnetic field. The latter is related 
to the electromagnetic tensor F^ v by — A Uyfl — A^^. 
Unless otherwise stated, we will use geometrical units, 
G = c = l. 
The Hamiltonian equations are given as 



d% d7r u 



dX 



dU 



(2) 



where A = r/m is the affme parameter, r denotes the 
proper time, and p^ is the standard kinematical four- 
momentum for which the first equation reads p^ = tt^ — 
qA* 1 . 

In the case of stationary and axially-symmetric sys- 
tems, we identify two constants of motion, namely, the 
energy E and angular momentum L. From the second 
Hamiltonian equation (j2j) we obtain 



7r t =pt + qA t = 
7T0=P0 + qA^ 



-E 
L. 



(3) 
(4) 



The trajectory is specified by the integrals of motion 
E, L, and the initial values r(0), 0(0) and u r (0). The 
initial u°(0) can be calculated from the normalization 
condition, g^ v u^u v = — 1 (we always choose the non- 
negative root). 

The effective potential can be derived in the following 
form: 



-J3 + a//? 2 - 4ct 7 
2a 



where 



a = -<r, 

l3 = 2[g^(L-qA 4> )-g tt qA t ], 
g^{L-~ q A,f-g tt ~eA\ 
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(5) 

(6) 
(7) 
(8) 



% = \g^{K fl -qA tl ){K v -qA v ), 



+2g t +qA t (L - qA^) - 1, 



where m and q are the rest mass and charge of the test 
particle, tt^ is the generalized (canonical) momentum, 



and where we introduce specific quantities L 



E 



and the specific charge q = j-. Local minima 
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Figure 3. Regular motion in the equatorial potential lobe in the fully integrable system of charged test particle (E = 0.99, L = 5M, 
q = 10 4 , r(0) = 32.02 M, 0(0) = 1.54) in the pure Kerr-Newman spacetime (Q = 3 x 10 -5 , a = 0.5 M) endowed with the fourth Carter 
constant of motion C. Long diagonals parallel to the LOI are general hallmark of regularity in the RPs. 
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Figure 4. Regular off-equatorial motion of a charged test particle (E = 1.77, L = 
Kerr background (a = 0.5 M) enriched with the Wald test field (qBo = 2M _1 , qQ ■ 
integrable systems are preserved, though the pattern is more complicated. 



time 

5M, r(0) = 2.9 M, 0(0) = 0.856 and u r (0) = 0) on the 
= 2). The diagonal structures typical for trajectories in 



of V e fi(r,0) reflect the location of stable orbits of test 
particles. Off-equatorial potential minima were identi- 
fied and vario us types of potenti al lobes were discussed 
elsewhere (see Kovar ^et al.ll2010l ). We can express the 
effective potential (j5j) as a function of r and u r and 
use it to determine the boundaries of allowed regions in 
Poincare surfaces of section for a given value of 00 

We employ the Kerr metr ic in standa rd Boyer- 

Lindquist coordinates t, r, 0, <j) (Misne r^t al.lll973j ): 



ds 2 



2 n 



a sin c 



(9) 



sin 



(r 2 +a 2 )d0-adt 



A 



£d<9 2 



1 We have employed the method of effective potential to study 
the stabilit y of the motion, however, we note that the force formal- 
ism (Abramowicz et al. 1995; Kovaf Sz Stuchlik 2007) can serve as 
a very efficient alternative tool. In particular, the off-equatorial 
motion of charg ed particles ca n be examined via the procedure de- 
scribed in Kovaf et al. (2010). This allows us to localize minima 
of the effective potential around which the stable orbits occur. 



where a stands for the spin parameter (the specific an- 
gular momentum, < a < M), A = r 2 — 2Mr + a 2 , and 
E = r 2 +a 2 sin 2 9. It is sufficient to consider positive val- 
ues of a without loss of generality (the cases of prograde 
and retrograde motion are distinguished by the sign of 
the particle charge and the orientation of the magnetic 
field). By setting a = the metric reduces to the static 
one of the Schwarzschild spacetime. 

At this point, a note is worth on the adopted compu- 
tational scheme which we have employed to study the 
trajectories and to detect the chaotical behavior. In or- 
der to reach reliable results, we coded several approaches 
and we checked their stability and precision. We employ 
the multi-step Adams-Bashforth-Moulton solver to de- 
termine the phase-space trajectory by numerical integra- 
tion of eqs. (j2j). In some cases, when a higher precision 
is demanded, we use the 7-8th order Dormand-Prince 
method that belongs to the family of explicit Runge- 
Kutta solvers with adaptive stepsize. This method im- 
proves the accuracy significantly, as can be verified by 
checking the conservation of the integrals of motion along 
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Figure 5. A transitional state between the regular and chaotic regimes of motion of a highly charged test particle which only differs from 
the previous case by increasing the energy to E = 1.796. The diagonal lines in the RP are partially disrupted, indicating the onset of chaos. 
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Figure 6. The chaotic motion of a highly charged test particle which only differs from the previous case by increasing the energy to 
E = 1.7975. The diagonal lines in the RP are now disrupted and complex large-scale structures appear which are a characteristic indication 
of deterministic chaos. 

tories are integrable, and so the chaos can set in only 
when perturbations of the background gravitational field 
are introduced or additional electromagnetic interaction 
with fields of external sources are allowed. This is also 
where recurrence analysis can be helpful. 

Methods of phase space recurrences have been success- 
fully applied to a wide range of various empirical data, 
not only in physics but also related to physiology, ge- 
ology, finances and other fields. Recurrence plots are 
especially suitable for the investigation of rather short 
and nonstationary data. On the other hand, the method 
of recurrence analysis has not yet been widely applied to 
study the dynamical properties of motion in relativistic 
systems. We thus briefly summarize this approach for 
our context. 

Besides more traditional methods of the numerical 
analysis of dynamical systems, such as a visual survey 
of Poincare surfac es of section or the evaluation of the 
Lyapunov spectra (|Skokosll2QlQh , the recurrence analysis 
is a rather novel technique, based on the analysis of re- 
currences of the system into the vicinity of its previous 
states. 

Recurrence Plots (RP) are introduced as a tool of visu- 
alizing the recurrences of a trajectory in the phase space 
(Eck mann et al.| [l987). The method is based on exam- 



the trajectory. However, the improved accuracy comes 
at the expense of computational time, as the adopted 
Dormand-Prince scheme is more computationally de- 
manding than the Adams-Bashforth-Moulton solver. 

Furthermore, it is well-know that, when dealing with 
Hamiltonian systems, the most appropriate solvers are 
those which res pect the symple ctic nature of Hamiltonian 
dynamics (e.g. lYoshidal fl993). We therefore employed 
also the implicit Gauss-Legendre Runge-Kutta (GLRK) 
method, which is a symplectic scheme. Indeed, we con- 
firm that this code provides the most reliable results, es- 
pecially in the case of long-term integration. The differ- 
ence in the accuracy between GLRK and non-symplectic 
solvers reaches several orders of magnitude and it is gen- 
erally more apparent in the case of chaotic trajectories, 
as expected. However, the cost in terms of the computa- 
tional time is also non- negligible, and so we only use the 
GLRK method to achieve very accurate long-time deter- 
mination of the trajectory in several exemplary runs. 

3. RECURRENCE ANALYSIS 

The Kerr metric is well-known and the analysis of test 
particle motion in this spac etime was carried out in many 
papers (Mis ner et al.l Il973h . Among important features 
of the Kerr metric is the fact that the particle trajec- 
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Distance to the nearest recurrence point [arbitrary unitsl 

Figure 7. Color map of the distance in the phase space to the nearest recurrence point. This example concerns a charged particle 
trajectory near the Kerr black hole in the asymptotically uniform magnetic field. Left: the case of regular motion with the energy of 
E = 1.77. Right: the case of chaotic motion with E = 1.7975. In the latter case more complex structure appear in the recurrence plot. The 
common parameters of both panels are L = 5M, a = 0.5 M, qBo = 2M _1 and qQ = 2 with the initial condition r(0) = 2.9 M, 6(0) = 0.856 
and u r (0) = 0. 

ination of the binary values that are constructed from 
the trajectory x(t). Results of the orbit analysis can be 
quantified statistically in terms of the Recurrence Quan- 
tification Analysis (RQA)0 

The RP construction is straightforward regardless of 
the dimension of the phase space. We only need to evalu- 
ate the binary values of the recurrence matrix , which 
can be formally expressed as follows: 

R i ,-(e) = e(e-||aj(i)-aj(j)ll), i,j = l,...,N, (10) 

where £ is a pre-defined threshold parameter, the 
Heaviside step function, and N specifies the sampling 
frequency. The sampling frequency is applied to the time 
segment of the trajectory x(t) under examination. There 
is, however, no unique prescription for the appropriate 
definition of the phase space norm || . || in eq. ([TO]) . We 
can consider a purely abstract vector space and apply 
one of the elementary norms L 1 , L 2 (Euclidean norm) or 
L°° (maximum norm). Some aspects of the appropriate 
choice of the norm are deferred to Appendix [A] 

Finally, we need to specify the value of the threshold 
parameter e. To th is end we follow the suggestion of 
M arwan et al.l ([2007, sec. 3.2.2) and relate e to the stan- 
dard mean deviation, a, of the given data set. Setting 
e = ka is advantageous because the proportionality con- 
stant fc, once adjusted to obtain a properly filled Recur- 
rence Plot, remains valid (with only minor adjustments) 
for all data sets of other trajectories. We therefore nor- 
malize the time series of each coordinate separately to 
zero mean and a = 1. 

The binary valued matrix R^ represents the RP which 
we get by assigning a black dot where R^ = 1 and leav 



(LOI). 

Recurrence Plots contain wealth of informati on about 
the dynamics of the system (|Thiel et al.ll2004gj ). Differ- 
ent pieces of knowle dge are encoded in large-scale and 
short-scale patterns (Mar wan et al.l [20071 sec. 3.2.3). To 
decide if a particular trajectory is a regular or a chaotic 
one, the determining factor is the presence of diagonal 
structures in the RP. Diagonal lines in the RP reflect the 
time segments of phase space trajectory during which the 
system evolution proceeds in a regular way. It captures 
the epoch when the trajectory proceeds almost parallel 
to its previous segment, i.e. within the £-tube around 
that segment. Hence, integrable systems exhibit them- 
selves by diagonally oriented structures in their RP. On 
the other hand, if the motion is chaotic the diagonal lines 
disappear and the diagonal features become shorter, as 
the trajectories tend to diverge quickly. As a result, more 
complicated structures appear in the RP. 

We stress that the interpretation of the RP is primar- 
i ly in tuitive. We refer to the review paper iMarwan et all 
(2007, sec. 3.2.3) where the patterns appearing in the 
RP and their relation to the current dynamic regime are 
analyzed in detail. Although basic conclusions may be 
inferred in general (e.g. distinction between regular ver- 
sus chaotic regime) the fine structure of the RP depends 
heavily on the properties of a given dynamic system. In 
order to gain more insight into the way in which vari- 
ous dynamical regimes manifest themselves in our system 
we typically present RPs accompanied by corresponding 
Poincare surface of section throughout this paper. 

Visual behavior of RP and it's complexity is quantita- 
tively reflected in RQA. The RQA evaluates statistical 
ing~a white dot where = 0. Both axes represent a characteristics of the recurrence matrix First of all, 



time segments over which the data set (the phase space 
vector) is being examined. RP is thus symmetric; the 
main diagonal is always occupied by the line of identity 



we define the recurrence rate (RR) as a density of points 
in RP, 



2 We use the CRP ToolBox (IMarwan et al.l [20071 . p. 321) in 
Mat lab (R2009b) to construct RPs and to evaluate RQA measures. 



RR(e) 



1 N 
N2 E R *.i( £ )- 



(11) 
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Figure 8. Graphs of different RQA measures based on the diagonal lines in the RP as a function of specific energy E. In each panel, 
400 trajectories were analyzed in a given energetic range. All measures exhibit the evident change of their behavior at E ^ 1.7954 which 
we interpret as an onset of chaos. Other parameters remain fixed at following values: L = 5M, qBo = 2, r(0) = 2.9 M, 6(0) = 0.856, 
u r (0) = 0, a = 0.5 M and qQ = 2. 



Now we can turn our attention to diagonal segments in 
RP. Their length draws distinction between regularity 
and chaos. The histogram P(e, I) records the number of 
diagonal lines of length I. It is formally given as follows: 



and the corresponding divergence (DIV) is defined as in- 
verse of the length of the longest diagonal line L m ax, 



DIV = 



1 



N 



(15) 



P(e, I) = J2 (! " Ri-i,i-i(e))(l - *>i+ij+i(e)) (12) 

l-l 



This histogram defines the determinism factor (DET), 
defined as a fraction of recurrence points, which form 
the diagonal lines of length at least l m [ n to all recurrence 
points, 



DET = 



(13) 



The average length of diagonal lines L (where only lines 
of length at least l m { n count) is 



DIV is in its very nature closely related to the diver- 
gent featur es of the phase space trajectory, and so it was 
originally (Eckm ann et al.lll987[ ) claimed to be directly 
related to the largest positive Lyapunov characteristic 
exponent A max . On the other hand, theoretical consid- 
erations justify the use of DIV as an estimator only for 
the lower limit of the sum of the positive Lyapunov ex- 
ponents (Ma r wan et al.ll2Q07l sec. 3.6). Nevertheless, a 
strong correlation be tween DIV and A max arises in nu- 
merical experiments (iTrulla et al .11 19961 ). 

The quantification measure ENTR is defined as the 
Shannon entropy of the probability p(e,l) = P{s,l)/Ni 
of finding a diagonal line of length I in the Recurrence 
Plot, 



ENTR = - ]T p(e,Z)]np(e,Z), 



(16) 



lP(e,l) 



l = ln 



(14) 



where Ni is a total number of diagonal lines: Ni(e) 
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Figure 9. Shannon entropy of probability distribution of diagonal lines lengths ENTR and three RQA measures based on the vertical 
lines of RPs as a function of specific energy E (details in the text). To some surprise, the vertical measures also react dramatically to the 
onset of chaos at E = 1.7954. 

c ) can 



Analogous statistics may be performed for vertical as 
well as the horizontal segments (RP is symmetric with 
respect to the main diagonal) . These segments are gener- 
ally connected with periods in which the system evolves 
during its laminar state. To this end, the histogram 
P(s, v) records the number of vertical lines of length v 
and it can be constructed as follows: 



N 



P(e,v)= J2(l-R id (e))(l 



v-l 
k=0 

(17) 

In analogy with the diagonal statistics histogram, 
P(e, v) is used to define the vertical RQA measures. 
Laminarity (LAM) is defined as a fraction of recurrence 
points that form vertical lines of length at least v m i n to 
all recurrence points, 



LAM ; 



(18) 



The trapping time (TT) is an average length of vertical 
lines, 



TT 



T. V Z7 m .P(e,v) 



(19) 



Finally, the length of the longest vertical line (V^ 
also be of some interest. 

RQA measures the crucial dependence of RP on the 
value of the threshold parameter, e, which must be ad- 
justed appropriately to a given data set. This lack of 
invar iance is a drawb ack of both RPs a nd RQA. Never- 
theless, it was shown ([Thiel et al.ll2QQ4bh that stable esti- 
mates of various dynamical invariants, such as the second 
order Renyi entropy and the correlation dimension, can 
be inferred if e is kept within a reasonable range. Since 
we shall use the standard RQA measures to compare the 
dynamics between test particles with different initial con- 
ditions, we have to eliminate the numerical effect of vari- 
ances in the range of coordinate values spanned by these 
trajectories. We achieve this by fixing the value of e. 

After the brief set of preliminaries we are now prepared 
to proceed to the intended application of the recurrence 
analysis in the next section. 

4. KERR BLACK HOLE IN UNIFORM MAGNETIC FIELD 

Large-scale magnetic fields are known to be present in 
cosmic conditions. They can exist around black holes, 
which do not support their own magnetic field but may 
be embedded in fields of distant sources. In the case of 
neutron stars, dipole-type magnetic fields of very high 
strength often arise. We concentrate on black holes in 
this section and defer the case of a magnetic star to sec. 
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Figure 10. Comparison of purely off-equatorial trajectories in spacetimes differing by the spin parameter a (left panels: a = 0.3M; 
middle: a = 0.6M; right: a = M) which is linearly linked to the energy E (left panels: E = 1.56; middle: E = 1.9, and E = 2.35 in 
the right panels). Other parameters remain fixed: L = 5M, M _1 , 6(0) = ^section = 0.856, qBo = 2M _1 and qQ = 2. RPs are taken for 
trajectories with r(0) = 2.9 M and u r (0) = 0. Although the structures in the Recurrence Plots clearly differ from each other, all of them 
represent diagonally oriented patterns that are characteristic of regular motion. The regularity of the motion is confirmed by surfaces of 
section in the bottom panels, where several trajectories (for each value of spin and energy) are presented. Different colours (grey scaled in 
the printed version of the paper) are used to distinguish the orbits originating from different initial conditions. 



El . , 

By employing the uniform test field solution (|Waldl 
119741 ) we incorporate a weak large-scale magnetic field 
near a rotating black hole. The vector potential can be 
expressed in terms of Kerr metric coefficients ([9]) as fol- 
lows, 



A t -- 
A 6 -- 



2ag t t) - 
- 2ag t(f) ) 



\Qgu — \( 



(20) 
(21) 



where Bo is magnetic intensity and Q stands for the test 
charge on the background of Kerr metric. The terms 
containing Q can be identified with the components of 
the vector potential of Kerr-Newman solution (although 
the test charge does not enter the metric itself). An 
example of an integrated trajectory is shown in Fig. [TJ 
IWaldl (j!974f ) has shown that the black hole selectively 
accretes charges from its vicinity, until it becomes itself 
charged to the equilibrium value 



Qw — 2Boa. 



(22) 



We remark that the particle charge q appears always 
as a product with Q or Bo in the formula (j5]) for the 
effective potential, as well as in equations of motion (|2j). 
Therefore, the simultaneous alteration of g, Q and Bq 



values, preserving the products qQ and qBo, does not 
affect the particle dynamics. If we further assume that 
Q — Qw = 2Boa is maintained, we only need to specify 
the value of qBo to uniquely determine a particular tra- 
jectory. However, since we do not restrict ourselves to 
the case Q = Qw we decide to always explicitely specify 
the values of qQ and qBo. 

4.1. Motion within the potential lobes 

Our previous analysis (jKovaf et al.l l2010f ) concludes 
that the off-equatorial bound orbits are allowed only for 
test particles obeying simultaneously the two conditions, 
sgn(aL) = 1 and sgn(g) = sgn(a£?o). Three distinct 
types were found (see Fig. [2j) and we examined the dy- 
namics of test particles in all of these types. The results 
for different types are comparable, and so we present 
here only the analysis of one of them, namely the type 
la. While raising the energy level from the local minima 
of the symmetric halo orbits, we observe that the off- 
equatorial lobes grow and eventually merge with each 
other once the energy of the saddle point in the equato- 
rial plane is reached. 

The size of the lobes is controlled by the specific en- 
ergy E. Employing the Poincare surfaces of section and 
the Recurrence Plots we investigate how the regime of the 
particle motion changes with E. This parameter appears 
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Figure 11. Comparison of trajectories of particles launched from the equatorial plane with different spin values. Also in this case we 
have to link linearly the value of spin a with E in order to maintain the existence of the potential lobe. In left panels we set a = 0.5M, 
E = 1.795, in middle panels a = 0.6M, E = 1.92 and in right panels a = M, E = 2.42. For all three cases we show surfaces of section 
of several trajectories differing in initial values r(0) and u r (0). Recurrence Plots are taken for trajectories with r(0) = 2.15M, it r (0) = 0. 
Other parameters remain fixed: L = 5M, qBo = 2M _1 , 6(0) = ^section = \ , qQ = 2. 



as a suitable control parameter producing a sequence of 
bound trajectories while all other parameters (and the 
initial position) remain fixed. For the sake of compar- 
ison we first present the case of a fully integrable sys- 
tem of a charged test particle on the pure Kerr- Newman 
spacetime (Fig. [3]). The motion occurs in the potential 
lobe around the local minimum in the equatorial plane; 
there are no halo orbits above th e horizon in this case 
(jKovaf et al]l2010t Ide Felicd[l979h . 

Figure 2] shows regular motion occurring in the off- 
equatorial lobe on the Kerr background with Wald test 
field. Increasing the energy level (while keeping all other 
parameters fixed) above the value in the equatorial sad- 
dle point results in a transitional regime depicted in Fig. 
El The onset of chaotic features does not occur as a di- 
rect consequence of the lobe merging. We rather observe 
that the orbit bound in merged (cross-equatorial) lobe 
remains regular until the particle notices the possibility 
of crossing the equatorial plane which happens when its 
energy is increased sufficiently. Once the motion becomes 
cross-equatorial, chaotic features appear. By increasing 
the energy even more we approach the critical value when 
the lobe opens and allows the particles to fall onto the 
horizon. For energies slightly below this limit we detect 
a fully chaotic regime of motion (Fig. [6]). 

Figure [71 shows an alternative representation of the re- 



currence plots, where the distance to the nearest neigh- 
boring point along the phase space trajectory is encoded 
by different colors. Again, by comparing the two pan- 
els one can clearly recognize how the diagonal structures 
disintegrate into scattered points as the chaos sets in. 

From the survey of Poincare surfaces of section and the 
Recurrence Plots we may only conclude that the transi- 
tion from the regular to chaotic regime occurs somewhere 
close to the value E = 1.796. In order to localize this 
transition more precisely we evaluate RQA measures for 
400 trajectories with the energy spread equidistantly over 
the interval E e (1.7948, 1.7968). In figs. flgj) and ®, we 
observe a sudden change of the behavior of the statistical 
measure at E « 1.7954, reflecting a dramatic change of 
the particle dynamics. 

Moreover, we know that the divergence DIV is related 
to the Lyapunov exponents, and in Fig. [8] we observe 
that it suddenly rises at E « 1.7954, meaning that the 
trajectories become more divergent when this energy is 
reached. All of these indications combined lead to the 
conclusion that this energy level represents a critical 
value at which a transition from regular to chaotic regime 
occurs. 

We conclude that the energy of the particle E acts 
as a governing factor determining the dynamic regime 
of motion. Our survey across various initial conditions 
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Figure 12. Profiles of the effective potential V e R for L = 5M, a = 0.5 M, qQ = 1.03, qBo = 0.1M -1 , taken in the equatorial plane and 
in the asymptotic region z = rcosO — > oo (as indicated with the corresponding curves). The valley crosses the equatorial plane, almost 
unaffected in its bottom parts, although the behavior of the potential near the symmetry axis is quite different in the equatorial plane 
where it approaches the horizon of the black hole (vertical dotted line at r = r+(a)). The horizontal dotted line at unity measures the rest 
energy of the particle m. 



has shown that motion in potential wells of the type la 
in a Wald test field is generally regular. Chaos appears 
well after the merging point of the lobes and close to the 
critical breaking energy. We have verified that all queried 
RQA measures, i.e. not only those based on diagonal lines 
in RP, react to the transition from the regular to chaotic 
regime of motion, allowing us to localize precisely the 
transition. This was fully confirmed also in the other 
two types (classes lb and Ic of our typology in Fig. [2]). 

4.2. The effect of spin on the chaoticness of motion 

Much attention has been recently focused towards the 
problem of determining the black hole spin from the 
properties of motion of surrou nding matter (Na ravanl 
120051 : iReynolds fe Nowall2QQ3f ). The astrophysical mo- 
tivation to address these issues arises from the fact that 
cosmic black holes are fully described by three parame- 
ters - mass, electric charge and spin. While the methods 
of mass determination have been widely discussed (e.g. 
Casaresl 120071 : iVestergaardl 120101 : ICzernv fc Nikolaiukl 



20l0h . the electric charge is considered to be negligible 



because of rapid neutralization of black holes via selec- 
tive accretion. However, determining the spin is a much 
more challenging task: the spin is important, but its in- 
fluence is apparent o nly very near the black hole horizon 
([Murphy et al.ll2009f ). 

One can raise a question of whether the value of spin 
parameter a of the Kerr black hole affects the dynamical 
regime of motion in the immediate neighborhood of the 
black hole. In other words, we ask if the spin parameter 
a triggers or diminishes the chaoticness of the system. 
Answering this question is not straightforward because 
by altering the spin across an interval of values (a 2 < 
M 2 ) we inevitably have to change some other variables 
of the system; otherwise the different cases could not 
be directly compared. Moreover, the location and the 
very existence of the potential lobes is not automatically 
ensured over the whole range of spin because of strong 
V e fi(a) dependence. 

We found that in order to keep the off-equatorial lobe 



at the initial position and the original size we need to 
increase the energy E, roughly proportionally to the in- 
crement of a. The effect of increasing a exhibits itself 
by lifting the hyperplane of effective potential. To com- 
pensate for this effect we have to elevate the i^-plane at 
which we cut the potential, so that we obtain (roughly) 
the original closed contour (the potential lobe) inside of 
which the motion is confined. It turns out that this can 
be achieved by linking both quantities linearly. 

First we compare the dynamics in the off-equatorial 
lobe in the range G (0.3, 1) (for a < 0.3M the topol- 
ogy of the effective potential changes) to which we lin- 
early relate the energy range E G (1.56,2.35) (whilst 
other parameters are kept fixed as follows: L = 5M, 
qB = 2M"\ r(0) = 2.9 M, 0(0) = 0.856, u r (0) = 0, 
qQ = 2). By inspecting the Poincare surfaces of sec- 
tion and performing the recurrence analysis for a large 
number of trajectories across the given range of a and 
E (exemplary cases presented in Fig. [TO]) we come to 
the conclusion that there is no overall trend that could 
suggest that a is the unique driving agent affecting the 
regime of motion. All trajectories in our survey exhibit a 
regular behavior, which is also in agreement with the pre- 
vious conclusion that the motion in off-equatorial poten- 
tial lobes associated with the Wald test field is generally 
regular. 

We also examined the dynamics of test particles 
launched from the equatorial plane whose trajectories 
occupy the potential lobe extending symmetrically above 
and below the equatorial plane. A given lobe maintains 
its size for spin values G (0.5, 1) and the related in- 
terval of energy E G (1.795,2.42). A survey across the 
given range of spin (energy) values reveals for this class of 
trajectories both chaotic and regular regimes. In Fig. [TT] 
we observe that for the lowest inspected spin, a = 0.5M 
(E = 1.795), the regular motion dominates, although is- 
lands of chaotic behavior are also present. Increasing the 
spin (energy) we observe that regular trajectories grad- 
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Figure 13. An exemplary trajectory (E = 1.058, L = 5M) is launched from the equatorial plane 0(0) = ^ with u r {0) = 0. Parameters of 

the background are a = 0.5M, qBo = 0.1M -1 , qQ = 1.03. In the upper left panel we observe that setting r(0) = 8 AM results in oscillations 
around the equatorial plane while launching it at r(0) = 8.7M makes it escape. In the upper middle panel we examine the trajectory of the 

escaping particle in terms of the rescaled radial coordinate r* = r r+ . In the case of oscillating trajectories two distinct modes of motion 
are observed (upper right panel). The first particle (r(0) = 11. 5M) shows a complex "ribbon-like" trajectory; the other one (r(0) = 8 AM) 
fills uniformly the given portion of the potential valley. The Recurrence Plots are also shown (bottom panels). We observe a highly ordered 
regular pattern for the particle with r(0) = 8AM (left panel), a more complicated diagonal pattern of the ribbon-like trajectory (launched 
at r(0) = 11. 5M, middle panel), and a disrupted diagonal pattern of the transitional trajectory (r(0) = 11. 4M, right panel). 



ually diminish. For a = 0.6M (E = 1.92) some regular 
orbits still appear, but they are already dominated by 
chaotic trajectories. For higher spins the traces of regu- 
lar motion further diminish. We present the extreme case 
a = M (E = 2.42) in Fig. [TT]to illustrate this apparent 
chaotic takeover. 

We conclude that for the class of orbits originating 
in the equatorial plane, spin a could possibly act as a 
destabilization factor which triggers the chaotic motion 
if enhanced sufficiently. However since the energy E is 
increased simultaneously it is not possible to attribute 
the observed dependence to the spin itself. On the other 
hand we have already seen in the above-given discussion 
(sec. 14. ip that energy E may itself act as a key factor 
determining the dynamic regime of motion. Thus we 
suggest to attribute the observed triggering of the chaos 
to the increase of energy E rather then spin a. 

We remark that a similar problem concerning the spin 
dependen ce of the motion chaot i cness was addressed re- 
cently by Takah ashi fc Koyamal ((2009). Authors of the 
quoted paper employ a dipole magnetic test field upon 
a Kerr background and perform a study of test particle 
trajectories, concluding that increasing the value of the 
spin parameter stabilizes the motion in a given setup. 
Unlike our case, the topology of the effective potential 
in their scenario allows the chosen potential lobe to be 



maintained at a given location and roughly the same size 
(but not the same depth) even if a varies while E is kept 
constant. However increasing the spin allows the authors 
to set gradually lower and lower energies, for which more 
regular trajectories are found. This is not at all sur- 
prising in perspective of our results, where the energy 
E proved to play a key role in determining the stability 
of motion. Although a direct comparison of presented 
surfaces of section differing only in a value may suggest 
the spin dependence, there is no clear and unambiguous 
correlation. We thus suggest attributing the observed 
dependence primarily to the level of energy, which acted 
as a motion destabilizer also in our setup. In fact, even 
if a real trend with the spin is present, it is hard to dis- 
entangle it from a simultaneous change of E. 

The question of the dependence of the dynamics upon 
the other parameters of the system (L, qBo, qQ) was also 
addressed during the analysis. It appeared that above 
mentioned difficulties accompanying the analysis of the 
spin dependence became even more serious in this case. 
Namely, neither it was possible to maintain the given po- 
tential lobe for a reasonable range of values of selected 
parameter nor we were able to fix this problem by bind- 
ing this parameter in some simple manner to some other 
parameter (e.g. energy E). In other words, none of these 
parameters itself may be regarded as a trigger for chaos. 
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Figure 14. The Poincare surface of section of several trajectories (E = 1.058, L = 5M, qB = 0.1M -1 , qQ = 1.03, a = 0.5M, 6 = ^ ) 
launched from the equatorial plane with various values of r(0) and u r (0) = 0. Grey colour indicates the escape corridor which lets the 
particles escape from the equatorial plane (see details in the text). 
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Figure 15. Selected types of the effective potential behavior in the 
with rotating dipole magnetic field. The inner bold line signifies the 

4.3. Motion in potential valleys 

Besides off-equatorial lobes, the potential may form 
another remarkable structure - an endless potential val- 
ley of almost constant depth which runs parallely to the 
symmetry axis (Fig. [T2]) . The poloidal orientation and 
asymptotical form of the valley are due to the fact that 
the Wald test field does not vanish at spatial infinity and 
it approaches the uniform magnetic field parallel to the 
symmetry axis. 

The existence of such a potential corridor suggests that 
test particles with a particular range of parameter values 
could escape from the equatorial plane to large distance. 
We observe that a test particle in the potential valley can 
keep oscillating around the equatorial plane or it may es- 
cape from this plane completely, depending on its initial 
position in the phase space (see the upper left panel of 
Fig. IT3|) . We can examine the motion in the asymptotic 
region by rescaling the radial coordinate (upper middle 
panel of Fig. [13]). To achieve this, we use r* = r ~ r+ , 
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r sine [M] 

vicinity of off-equatorial halo orbits above the surface of magnetic star 
surface of the star at r = 4M. The outer line is the light surface. 



where r + = M + \J M 2 — a 2 is the position of the outer 
horizon of the black hole. 

The normalization condition allows one to express u e 
as a function of other phase space variables and param- 
eters of the system. Intuition suggests that, considering 
particles launched from the equatorial plane (6(0) = §), 
the initial value u e (0) is a governing parameter which 
decides whether the particle remains oscillating around 
= 77 or leaves it once forever. We can draw the isolines 

of selected u e values in the (r, u r )-plane which we use 
as a surface of section (0 sec tion = f ) f° r the inspection 
of the test particle dynamics. Comparing acquired iso- 
lines for various values of u e with the empirically stated 
escape corridor of Fig. [14] leads to the conclusion that 
they never coincide perfectly, although the correlation is 
quite high. In other words there is no definite threshold 
value of u e (0) which would determine whether a selected 
combination of r(0), u r (0) (while other parameters are 
fixed) lets the particle launched at equatorial plane leave 
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Figure 16. Regular motion in the off-equatorial potential lobe at the energy level E = 0.8482. Parameters of the system are qA4 = 
-5.71576 M 2 , L = 0.87643 M, Q = 0.011485 M _1 . The left panel shows sections of several trajectories launched from 6(0) = Section = 
1.0492. One of them (r(0) = 5.02 M and u r (0) = 0) is visualized in the Recurrence Plot in the right panel. The motion is regular (RP 
remains diagonal), however, the density of recurrence points clearly grows during the analyzed period. 



or oscillate around the plane. 

In Fig. [14] we observe several qualitatively different 
types of possible particle dynamics. Next to the effec- 
tive potential contour we find closed and well-defined 
curves which represent the regular motion. Going fur- 
ther inside we notice that fragmented curves are present. 
Below them we find progressively more and more blurred 
patterns, which again signifies the onset of chaos. The 
inner parts of the potential lobe are occupied by the "es- 
cape corridor" where the particles can stream freely from 
the equatorial plane, as seen in Fig. [131 

The trajectories represented by the closed curves in the 
surface of section differ profoundly from the disconnected 
curve orbits. The difference can be seen also in the direct 
projection onto the poloidal plane (upper right panel of 
Fig. [13]). While the trajectories of the first type gradually 
fill each particular compact region of given section, the 
latter forms bundles which curl through the projection 
plane resembling ribbons that bound regions which are 
never reached by the particle. 

By employing the Recurrence Plots (bottom panels of 
Fig. [13]) we confirm that the dynamics differs significantly 
in these distinct modes of motion. Not surprisingly we 
obtain a typical regular pattern in the case of trajectory 
which forms a closed sharp curve in the surface of section 
(bottom left panel of Fig. [13]). A "ribbon-like" trajec- 
tory (fragmented curve in Fig. [14]) results in an ordered 
checkerboard pattern (bottom middle panel), which is 
know n to be typical for pe riodic and quasi-periodic sys- 
tems (jMarwan et all [2007). Finally, in the bottom right 
panel of Fig. [13] we observe that a blurred curve tra- 
jectory exhibits slight chaotic behavior in its RP. The 
diagonal structures are partially disrupted and we notice 
that diagonal lines become bent as they approach the 
line of identity. 

Authors of the recent paper 

(jLukes-Gerakopoulos et al.l [2010) observed patterns 
similar to our fragmented curves in Poincare sections of 
phase space occupied by geodesic trajectories around 
slightly perturbed Kerr spacetime described by Manko- 
Novikov metric. They assume that such fragments can 
be identified with the Birkhoff islands anticipated by the 



Poincare-Birkhoff theorem. A given chain of Birkhoff 
islands (i.e. the fragmented curve in the section) 
originates from the resonant torus of an unperturbed 
system whose orbits have characteristic frequencies in 
a given rational ratio. The most prominent chains of 
islands should belong to the resonances characterized 
by simple integer ratios ^, | etc. Authors of the cited 
paper conclude that the presence of the Birkhoff islands 
might be in principle detected (at least in the context of 
gravitational radiation emitted by extreme mass ratio 
inspiraling sources) which would prove the system to 
be perturbed. However, since the assumption of the 
weakness of the perturbation is generally not fulfilled in 
our case, we do not explicitely identify fragmented tori 
in the surfaces of section with the Birkhoff chains. 

The class of escaping trajecto ries in the K err back- 
ground was recently discussed by iPretil (|2010f ) . The au- 
thor suggests that a Wald electromagnetic field employed 
in this setup could serve as a charge separation mecha- 
nism for astrophysical black holes since the sign of the 
particle charge may determine whether a given particle 
escapes from the equatorial plane or becomes trapped in 
cross-equatorial confinement (or falls into the horizon). 

5. A MAGNETIC STAR 

We describe the gravitational field outside a compact 
star by the Schwarzschild metric, 



ds2 = -[ 1 -—) dt2 - 

+r 2 (d<9 2 + sin 2 Odcf) 2 ). 



i-^yv (23) 

r J 



The associated magnetic fiel d is modeled as a dipole ro- 
tating at angular velocity Vt (Sengupta 1995|): 



A* = —SlAth = 



8M 3 



A, 
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where 



1Z = 2M 2 + 2Mr + r 2 log [ 1 - ^ j . 



(24) 
(25) 

(26) 
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Figure 17. For an energy value of E = 0.8485, both off-equatorial lobes merge via the equatorial plane. The upper-left panel shows 
Poincare sections of two trajectories (#(0) = # sec tion = 1.0492, u r (0) = 0). A particle launched at r(0) = 4.8 M never crosses the equatorial 
plane and moves regularly. Setting r(0) = 5 M we observe a chaotic motion crossing the equatorial plane repeatedly. All particles launched 
with r(0), u r (0), corresponding to the inner parts of the potential curve, move in the same chaotic manner. The outskirts are occupied 
by regular trajectories. The upper-middle panel shows the transient trajectory (r(0) = 4.85 M), regular during the integration period of 
A = 10 5 . In the upper-right panel, the integration time is prolonged to A = 3 x 10 5 . Here, the onset of chaos is connected with the first 
passage through the equatorial plane. The RP of the regular trajectory with r(0) = 4.8 M is presented in the bottom-left panel; the RP 
on the right belongs to the chaotic trajectory with r(0) = 5 M. 



The r elated dipole moment M is given by (Ba kala et al.l 
[20101 ) 



M 



4,M 3 r 3 J 2 (r* 



2M) 1/2 B 



simplifies to the form 
3qMlin sin 2 



eff : 



6M(r* — M) + 3r* (r* - 2M) In (l - 2Mr* 



(27) 

where Bo is the magnetic field at the neutron star equa- 
tor, r* is the radius of the star surfaceQ 
We assume eqs. (I23l) - (l24l) to hold outside the star sur- 



8M 3 
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(28) 
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r sin 6 



8M 3 r 



face (r > r*) and inside the light cylinder (u^u^ < 0). 
We set r* = 4M as the inner radial boundary of the 
particle motion. As for the light cylinder, the mentioned 
condition results in a relation r 2 sin^ 2 ^ 2 = 1 — which 
implicitly specifies the outer boundary. The vector po- 
tential (|211) is valid inside the rigidly corotating magne- 
tospheric plasma, which we consider to be an excellent 
conductor, so that the force-free condition F^u v = 
holds for the plasma for which — (u*,0,0,u^) and 

4 = o. 

A general formula for the effective potential eq. (j5j) 

3 The existence of extreme ly compact stars with r* 3M 
is unlikely, but not excluded (Bahcall et al. 1989; Stuchh'k et al. 
2009). Mos t of the realistic eq uations of state imply a lower limit 
r* ~ 3.5M (Glendenning 1997). On the othe r hand, the model s of 
Q-stars do allow a lower limit of r* w 2.8M ( Miller et al. 1998h . 



5.1. 



Motion inside the potential lobes 

Our previous analysis (jKovaf et al.l l2010f ) revealed a 
number of distinct types of possible topological struc- 
tures of the effective potential. However it appears that 
the system is not as rich in its dynamical properties. 
The test particle trajectories share some similar features 
across different classes of the effective potential. There- 
fore, we only present surveys of particle dynamics in 
three exemp l ary types: la, Ha and IIIc (Fig. [I5j see 
iKovaf et al.l (j2010f ) for the complete review). 

Class la lobes grow with energy increasing. Once the 
level of the equatorial saddle point is reached, the lobes 
merge with each other across the equatorial plane. The 
single merged lobe eventually intersects the surface of the 
star if the energy level is increased sufficiently, letting the 
particle fall onto the surface. Lobes of Ha type also merge 
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Figure 18. For the energy level E = 0.857 we obtain a broad lobe which almost touches the surface of the star at r = 4 M . In the 
upper-left panel, we launch two particles with 0(0) = 1.0492, r(0) = 4.75 M. The particle to the left of the star starts with w r (0) = 
and moves chaotically, while the other one with it r (0) = 0.03 follows a perfectly regular trajectory. The upper right-panel shows these two 
types of trajectory appear in the surface of section plot. The bottom panels demonstrate the difference between the two types in terms of 
Recurrence Plots. 



via the equatorial plane but in contrast to the first type 
the merged lobe opens toward the light cylinder (beyond 
which the model becomes invalid). Lobes of the class IIIc 
first open via the off-equatorial saddle points, allowing 
the particles fall onto the star, before the lobes merge 
through the equatorial plane. 

Now we study the three selected types of the effective 
potential topology in more detail. We are primarily in- 
terested whether and how the dynamic regime changes 
across the given range of specific energy E. Especially, 
we shall address what happens with the dynamics when 
the particle acquires enough energy to cross the saddle 
point. 

The first survey (type la) begins at energy level E = 
0.8482, corresponding to the closed lobe. In Fig. [TBI we 
observe that the motion inside the lobe is stable. No 
chaotic properties are detected - neither in Poincare sur- 
faces of section nor in the Recurrence Plots. 

As the energy increases to E = 0.8485, the symmetri- 
cal lobes merge via the equatorial plane. By inspecting 
a number of trajectories in this case we find that chaos 
starts appearing at this point - those particles which no- 
tice the gate through the equatorial plane always fill the 
entire allowed region and they move chaotically. Never- 
theless, there are still such particles which move regularly 
in one of the two parts of merged lobe and never cross 



the equatorial plane. We can also find transient trajec- 
tories corresponding to regular motion lasting for some 
period of time in one part of the lobe, followed by chaotic 
motion over the entire lobe once the particle finds and 
encounters the passage across the equatorial plane. All 
of the mentioned cases are illustrated in Fig. [T71 

Increasing the energy further to E = 0.857, we ob- 
tain a broad potential lobe which almost touches the 
star surface. The situation changes from the previous 
case where the gate connecting the off-equatorial lobes 
was narrow. Now we do not find trajectories which oc- 
cupy only one part of the lobe, above (or below) the 
equatorial plane. Chaotic trajectories densely filling the 
entire lobe are typical for this setup. We also encounter 
perfectly regular trajectories forming ribbon-like struc- 
tures spanned between northern and southern borders of 
the lobe. In Poincare sections these appear as regular 
islands surrounded by a chaotic ocean (upper panels of 
Fig. [15]). We notice that the RP of this regular trajectory 
is extraordinarily simple and consists of almost perfect 
diagonal lines (bottom panels of Fig. [T8|) . Thus its dy- 
namic properties are close to those of a periodic system, 
which is in contrast with the neighboring fully chaotic 
orbits. 

The second type (class Ha) of the effective potential 
topology of off-equatorial lobes differs from the first one 
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Figure 19. Regular motion in an off-equatorial lobe of the third type. Parameters used: E - 
L = 6.25382 M,Q = 0.011485 M _1 . Particles are launched from latitude 0(0) = 6> sec tion = 



= 0.99579, 
1.0492. 



L = 6.25382 M, qM = 45.87368 M 



significantly as the lobes do not open towards the star 
when the energy is raised sufficiently. On the contrary, 
in this case we observe that the lobe's boundary touches 
the light cylinder first. 

The motion in the off-equatorial lobes proves to be reg- 
ular while the merging lobes bring chaos into play. Chaos 
becomes dominant for broader lobes, however, stable reg- 
ular orbits also persist. The results are similar to those 
of the first type of potential topology discussed above. 

The last analyzed topology of the lobes (class IIIc) 
differs profoundly from the preceding two cases, as can 
be seen in Fig. [TSl We find that stable motion dominates 
in this setup. This can be verified by comparison with 

Fig.m 

As we further increase the energy level, we obtain more 
complicated shapes of the equipotentials that allow the 
particle to fall on the star surface. On the other hand, 
opening the outflow gate energetically precedes the merg- 
ing point of both off-equatorial lobes. In Fig. [20j we 
discuss the motion governed by the largest possible lobe 
which almost touches the light cylinder. We observe that 
stable regular orbits are still possible for those particles 
that do not hit the passage. 

From the above-given discussion we conclude that the 
motion of charged test particles in the off-equatorial lobes 
allowed by the test field of the rotating magnetic dipole 
on the Schwarzschild background is largely regular. Once 
the off-equatorial lobes merge with each other, chaos may 
appear. Increasing the energy, the chaotic motion be- 
comes typical but, quite surprisingly, very stable orbits 
also exist under these circumstances. 

6. DISCUSSION AND CONCLUSIONS 

We studied the regular and chaotic motion of electri- 
cally charged particles near a magnetized rotating black 
hole or a compact star. We employed the method of re- 
currence analysis in the phase space, which allowed us 
to characterize the chaoticness of the system in a quanti- 
tative manner. Unlike the method of Poincare surfaces, 
the Recurrence Plots have not yet been widely used to 
study the chaotic systems in the regime of strong gravity. 

The main motivation for these investigations is the 
question of whether the matter around magnetized com- 
pact objects can exhibit chaotic motion, or if instead the 



system is typically regular. One of the main applications 
of our considerations concerns the putative envelopes of 
charged particles enshrouding the central body in a form 
of a fall-back corona, or plasma coronae extending above 
the accretion disk. While we concentrated on the spec- 
ifications of the RP method in circumstances of a rela- 
tivists system, the assumed model cannot be considered 
as any kind of a realistic scheme for a genuine corona. 
We simply imposed a large-scale ordered magnetic field 
acting on particles in a combination with strong gravity. 

Various aspects of charged particle motion were ad- 
dressed throughout this paper. First of all, we investi- 
gated the motion in off-equatorial lobes above the hori- 
zon of a rotating black hole (modeled by Kerr metric 
equipped with the Wald test field), as well as above the 
surface of a magnetic star (modeled by the Schwarzschild 
metric with the rotating dipolar magnetic field). In both 
cases we conclude that the motion of test particles is reg- 
ular, which was confirmed for a representative number of 
orbits across the wide range of parameters over all topo- 
logical types of off-equatorial potential wells. This result 
is somewhat unexpected because the off-equatorial orbits 
require a perturbation to be strong enough (in terms of 
strength of the electromagnetic field), so that it can bal- 
ance the vertical component of the gravitational force. 

Further, we investigated the response of the particle 
dynamics when the energy level E was raised gradually 
from the potential minimum to values allowing cross- 
equatorial motion. We examined various topological 
classes of the effective potential and came to the conclu- 
sion that the cross-equatorial orbits are typically chaotic, 
although very stable regular orbits may also persist for a 
certain intermediate en ergy range. The classical work of 
iHenon fc Heilesl (j!964f ) should be recalled in this context 
since it also identifies the energy as a trigger for chaotic 
motion in the analysed simple system. More recently 
the Henon -Heiles system wa s revis ited in the relativistic 
context by IVieira fc Letelierl (jl996h . 

We also addressed the question of spin dependence of 
the stability of motion for Kerr black hole in the Wald 
field. We noticed that this is a rather subtle problem. 
The effective potential is by itself sensitive to the spin 
value a - hence, we had to link the potential value 
roughly linearly with the energy E to maintain the po- 
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Figure 20. For the energy level = 0.9962 (other parameters as in Fig. ll9|) we obtain a large lobe that almost touches the light cylinder 
and allows the particle to fall onto the surface of the star via a narrow passage above and below the equatorial plane. The upper left panel 
shows two trajectories launched from 6(0) = 1.0497, u r (0) = 0. One of the particles (starting from r(0) = 61.5 M) follows an unstable 
path and eventually falls on the star surface. On the contrary, the other particle (starting from r(0) = 72.5 M) moves regularly and never 
escapes any given part of the lobe. The upper-right panel shows these two kinds of trajectory depicted in the Poincare surface of section. 
Chaotically dispersed points belong to the escaping trajectory. The bottom-left panel shows the Recurrence Plot of the escaping particle; 
this plot does not exhibit typical chaotic behavior, although the large-scale structures are present. The bottom-right panel presents the 
Recurrence Plot of stable motion. 



tential lobe at a given position. In other words, we did 
not find any clear and unique indication of the spin de- 
pendence of the motion chaoticness. Most trajectories 
exhibited regular behavior, which is also in agreement 
with the previous results indicating that motion in off- 
equatorial lobes is generally regular. On the other hand, 
in the case of the cross-equatorial motion we observed 
that, for higher spins, more chaotic features come into 
play when compared with the case of slow rotation. This 
trend might be also attributed to simultaneous adjust- 
ments of E. In other words, it appears impossible to 
give an unambiguous conclusion about the spin depen- 
dence of the particles dynamics. Instead, one has to deal 
with a complex, interrelated dependence. 

In the case of a Kerr black hole immersed in a large- 
scale magnetic field, we observed the effect of confine- 
ment of particles regularly oscillating around the equa- 
torial plane. Escape of particles from the plane is allowed 
for a given range of initial conditions since the equipo- 
tentials do not close; they form an endless axial "valley" 
instead. The escaping trajectories create a narrow, colli- 
mated structure parallel to the axis. 

Asking the question about the dynamical properties of 
motion of matter in the black hole coronae is of the- 
oretical interest by its own, but it is relevant also in 
view of precise spectroscopical and timing studies with 



the forthcoming X-ray satellites. However, one more 
word of caution is appropriate especially concerning our 
assumption of the ordered magnetic field. Although 
the large-scale magnetic fields are adequate to describe 
the initial background field of the magnetic star or a 
toroidal current in the accretion disk ( Koide et al. 2006; 
iBeckwith et~alll2008l : iRothstein fc Lov elace 2001 e.g.), a 
turbulent compon ent is known to develop quic kly in the 
accreted plasma (McKin nev fc Naravanl [2007). These 
will also perturb the motion of matter, and so the regular 
orbits may quickly disappear. This should translate to 
the transient emergence and subsequent disappearance 
of the periodic component in the observed signal. 
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APPENDIX 
CHOICE OF PREFERRED OBSERVERS 

When analyzing the dynamics in the general relativistic context the fundamental questio n arises wheth er the dis- 
tinction between chaotic and regular motion is coordinate dependent or not. To this end Motter (2003) infers the 
transformation law for Lyapunov exponents. He concludes that although the Lyapunov exponents themselves are not 
invariant they transform in such a way that positive Lyapunov exponents remain positive and vice versa. In other 
words the distinction between regular and chaotic dynamics may be drawn invariantly. 

In order to gi ve the notion of recurrence a rigorous and, at the same time, an intuitive sense, we can employ the 
3 + 1 formalism (jThorne fc MacDonald 1982) that is based on an appropriate selection of a family of spacetime- filling 
three-dimensional spacelike hypersurfaces (foliation) of constant time t. The timelike curves orthogonal (in a spacetime 
sense) to the hypersurfaces may be regarded as the world-lines of a family of fiducial observers (FIDO) who naturally 
parameterize their world line by proper time r (whose rate of change generally differs from that of t). FIDO identify 
each spatial hypersurface along his world line as a slice of simultaneity. The geometry of this spacetime slice is given 
by 3-metric 7^: 

lij = 9ij + UiUj, (Al) 

where ui stands for the spatial part of FIDO's four- velocity and g^ for the spatial part of the spacetime metric. Con- 
sidering also the time coordinate, 7^ can be regarded as a projector onto the three-dimensional spatial hypersurface. 

We select zero-angular momentum observers (ZAMO) (Ba rdeen et al .1 11 9 72) as an appropriate representation of the 
fiducial observers. Their ort ho normal tetrad is given by 

O (A2) 

(A3) 
(A4) 
(A5) 

where A = (r 2 + a 2 ) 2 

Projecting an arbitrary four- vector C M onto ZAMO's hypersurface of simultaneity directly results in the following 
3-dimensional quantities, 

SB c i = 1 ij Cj = (c\C\^C^) , (A6) 

3B Ci = >y??C j = (C r: C e: g^C*) ; (A7) 

here, ^ (u^ for ZAMO) and, consequently, 7 < ^#</></> 1. 

ZAMO tetrad components 3D C^ are given as 

3D C (i) =3D c = e 3D C . = (^C"\ JgfeC\ y^Cf) . (A8) 



The hypersurface components of the phase space constituents x l and 7Tj, as measured by ZAMO, are then: 



e (t) = 


u = — 


f d 
K dt 


e (r) = 


VA d 
p dr 




e(8) = 


1 d 
~p"d6 




e (4>) = 


P d 
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= 2aA~ 


~ x Mr. 





3D_(i) _ 




(A9) 

1 ^ (A10) 



The spatial 3-metric is 

ds 2 = 5 m) dx®dxW + 0(\x^\ 2 )dx^dx^\ (All) 

where represents the spatial distance from the origin of the tetrad, i.e. ZAMO's current location. ZAMO is not an 
inertial observer, which generally causes the first order corrections 0(\x^ |) to the Minkowskian metric g(i)(j) = V(i)(j)- 
But these do not enter the spatial part of the metric. Thus the 3-metric within the spatial hypersurface is a Euclidean 
one, with the deviations of second order in the distance from the spatial origin on ZAMO's world-line. 
We will use ZAMO's metric at distances up to the value of the threshold parameter e. The Euclidean metric 

according to eq. (|A11|) will be therefore justified if <C 1, where P stands for a constant (for a given ZAMO). P 

has the dimension of length and characterizes the curvature of the hypersurface. We suggest setting P = if -1 / 4 , 
where K = R^ v ^ R^ v ^ represents the Kretschmann scalar evaluated from the Riemann curvature tensor. In the case 
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of the Kerr black hole the Kretschmann scalar may be expressed in a surprisingly simple form (Henry 2000). While 

2 

constructing the Recurrence Plots, we check whether the condition pj < 1 remains satisfied. 

The above-mentioned adoption of preferred observers is needed in order to maintain an operational criterion of 
chaos and be able to formulate an explicit form of the equations for RQA measures in a curved spacetime. Here the 
notion of the phase space distance plays a role. In Kerr metric (or another axially symmetric stationary spacetime), 
Fiducial Observers (a.k.a. FIDOs) represent a natural selection of preferred observers. Obviously, this option is not 
unique, and so a detailed appearance of the recurrence plots is also ambiguous to certain extent. But not so the main 
conclusions that we infer regarding the chaoticness versus regularity of the system behavior, because this distinction 
can be eventually traced down to the exponential versus polynomial growth of the separation with the particle proper 
time along neighboring trajectories. 

We can deduce the kind of transformation between different families of observers that could affect our conclusions: 
these are transformations involving exponential dependencies on observer's phase-space position. For example trans- 
formation to accelerated frames and spacetime points in the vicinity of singularities may need a special consideration, 
as well as the investigation of highly dynamical spacetimes that are lacking symmetries. On the other hand, se- 
lecting LNRF to define ZAMOs in (weakly perturbed) Kerr metric outside the black hole horizon appears to be a 
well-substantiated choice. 

Similar arguments for the adoption of pre f erred observers on the basis of spacetime symmetries have been elaborated 
in greater detail by iKaras fc Vokrouhlickyl (pj)92) in the context of Ernst's magnetized black hole, which is another 
particularly simple (static) exact solution of Einstein-Maxwell equations exhibiting the onset of chaos as the magnetic 
field strength is increased. 
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